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We describe the results of the two methods we developed to calculate the stationary nonlinear 
solutions in one-dimensional plasmonic slot waveguides made of a finite-thickness nonlinear dielectric 
core surrounded by metal regions. These two methods are described in detail in the preceding article 
[Walasik et al., submitted to PRA]. For symmetric waveguides, we provide the nonlinear dispersion 
curves obtained using the two methods and compare them. We describe the well known low-order 
modes and the higher-modes that were not described before. All the modes are classified into two 
families: modes with and without nodes. We also compare nonlinear modes with nodes with the 
linear modes in similar linear slot waveguides with a homogeneous core. We recover the symmetry 
breaking Hopf bifurcation of the first symmetric nonlinear mode toward an asymmetric mode and 
we show that one of the higher modes also exhibits a bifurcation. We study the behavior of the 
bifurcation of the fundamental mode as a function of the permittivities of a metal and a nonlinear 
core. We demonstrate that the bifurcation can be obtained at low power levels in structures with 
optimized parameters. Moreover, we provide the dispersion curves for asymmetric nonlinear slot 
waveguides. Finally, we give results concerning the stability of the fundamental symmetric mode 
and the asymmetric mode that bifurcates from it using both theoretical argument and numerical 
propagation simulations from two different full-vector methods. We investigate also the stability 
properties of the first antisymmetric mode using our two numerical propagation methods. 

PACS numbers: 42.65.Wi, 42.65.Tg, 42.65.Hw, 73.20.Mf 


I. INTRODUCTION 

Nonlinear slot waveguides (NSWs) are structures made 
of a finite-size nonlinear dielectric layer sandwiched be¬ 
tween two semi-infinite metal layers. They have been 
studied in Refs. mm where it was shown that they al¬ 
low for sub-wavelength confinement of light and phase 
matching for the second-harmonic generation. More re¬ 
cently, in the works of Rukhlenko et al. mm, analyti¬ 
cal formulas for the dispersion relations of these NSWs 
were presented for symmetric and antisymmetric nonlin¬ 
ear modes only. These dispersion relations were given us¬ 
ing integral equations that have to be solved numerically. 
The study of Davoyan et al. [5] showed, using the nu¬ 
merical shooting method to solve Maxwell’s equation in 
NSWs, that a symmetry breaking bifurcation that gener¬ 
ates an asymmetric mode from the fundamental symmet¬ 
ric mode occurs in NSWs. Such bifurcation phenomena 
are well known in fully dielectric nonlinear structures 0- 
[13]. Recently, higher order modes were also reported 
in NSWs [14]. Moreover, it was shown that plasmonic 
coupling and symmetry breaking phenomena can be ob¬ 
served in waveguides built of a linear dielectric core sand¬ 
wiched by nonlinear metals usms]. Nonlinear switching 
was predicted in NSW-based structures using numerical 
simulations m- 
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In the preceding article m, we describe two models we 
developed to study the stationary nonlinear solutions in 
NSWs where the nonlinear core of the focusing Kerr type 
was considered. The first model assumes that the non¬ 
linear term depends only on the transverse component of 
the electric field and that the nonlinear refractive index 
change is small compared to the linear part of the refrac¬ 
tive index. It allows us to describe analytically the field 
profiles in the whole waveguide. It also provides a closed 
analytical formula for the nonlinear dispersion relation. 
This first model is called Jacobi elliptic function based 
model (JEM). The second model takes into account the 
full dependency of the Kerr nonlinear term on all electric 
field components and no assumption is required on the 
amplitude of the nonlinear term. The disadvantage of 
this approach is the fact that the field profiles must be 
computed numerically. This second model is called the 
interface model (IM). 


This article is organized in the following way. In Sec.|TI| 
we describe the results obtained with our two models for 
symmetric NSWs. They include a mode classification 
taking into account the higher order modes we found 
previously El and a detailed study of the field profile 
transformation as a function of power. We also provide 
a permittivity contrast study that allows us to decrease 
by several orders of magnitude the bifurcation threshold 
at which the first asymmetric mode appears. In Sec. m 
we provide the results concerning asymmetric NSWs in 


which the mode degeneracy is lifted. Finally, in Sec. IV 
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using both theoretical arguments and numerical propa¬ 
gation simulations from two different full-vector methods 
we provide results on the stability of the main stationary 
solutions obtained in symmetric NSWs. 


II. RESULTS FOR SYMMETRIC WAVEGUIDES 

A. Dispersion relations, field profiles and mode 
classification 

In this section, the dispersion relations obtained for 
the symmetric NSW are presented. The field profiles 
corresponding to each of the dispersion curves are de¬ 
picted and allow us to classify the modes according to 
their symmetry and the number of nodes in the magnetic 
field profile. 

Figure presents dispersion relations for the sym¬ 
metric NSW obtained using the JEM and the IM. The 
parameters of the NSW studied are: €i = €3 = —90 
(gold), €1^2 = 3.46^, a 2 = 6.3 • 10“^^ m^/V^ (hydro¬ 
genated amorphous silicon) and d = 400 nm at a wave¬ 
length A = 1.55 jam. The geometry of the structure with 
its parameters is shown in Fig. The dispersion rela¬ 
tions present the dependence of the effective index of the 
mode /3 as a function of the power density in the waveg¬ 
uide core Pc which is calculated in the following way: 



FIG. 1. Dispersion diagrams l3{Pc) for the symmetric NSW 
obtained using (a) the JEM and (b) the IM. 
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FIG. 2. Geometry of the plasmonic NSW with the parameters 
of the structure. 


Pc= /%.dx, (1) 

Jo 

where Sz denotes the ^-component of the Poynting vector 
S = l/25?e(E X H*). 

We observe a very good qualitative agreement between 
the dispersion diagrams obtained using our two models. 
The number and the character of the dispersion curves 
is very similar in both cases. The qualitative agreement 
between the results of the two models confirms their va¬ 
lidity. Quantitatively speaking, the models agree in the 
range of low power densities (below 10^ W/m). Above 
this value we observe quantitative differences in the re¬ 
sults. The origin of the differences is explained by the 
assumptions made in the JEM (low nonlinearity, only 
Ex component of the electric field contributing to the 
Kerr nonlinear effect) as described in Ref. [18]. In the 
following, we will focus on the results obtained using the 
more accurate IM. 

The NSW supports numerous modes with various 
properties. Firstly, we will discuss the mode classifica¬ 
tion according to the symmetry of the mode. For the 
low power region, the NSW studied here supports two 
modes: a fundamental symmetric mode [blue curve la¬ 
beled SO in Fig. and blue field profiles in Fig. and a 
low-power antisymmetric mode [red curve labeled ANO in 
Fig. and red field profiles in Fig. [^. At Pc ~ 10^ W/m 
a symmetry breaking bifurcation occurs that gives birth 
to an asymmetric mode [5| [green curve labeled ASl in 
Fig. 0 and green field profiles in Fig. [^. This modes 
and this type of behavior are known in nonlinear waveg¬ 
uides [SHISlIIilllSl- The power density Pc of the modes 
SO, ANO, and ASl first increases with the increase of the 
effective index (3 and decreases for P ^ 4.75. 

Our models allow us to find new, higher order modes 
in NSWs. The higher order modes can be divided into 
two families: node-less modes and modes with nodes. 
Among the node-less modes we find higher order sym¬ 
metric modes (SI and SII) from which asymmetric modes 
bifurcate (AS2 and AS3, respectively). Their disper¬ 
sion curves are labeled with the name of the mode in 
Fig. I^and their field profiles are presented in Figs. and 
l^e), (f). Higher order node-less modes resemble a sin¬ 
gle soliton (SI and AS2) or two solitons (SII and AS3) 
propagating in the NSW core. 

All the dispersion curves of the asymmetric modes are 
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FIG. 3. Profiles of magnetic field Hy{x) for the symmetric SO 
mode (blue), antisymmetric ANO mode (red), and the hrst- 
order asymmetric ASl mode (green). The subplots present 
the transformation of the field prohles at the points corre¬ 
sponding to the vertical lines labeled a-e indicated in Fig.|^ 


FIG. 4. Prohles of magnetic held Hy{x) for the symmetric 
SI mode (blue) and the second-order asymmetric AS2 mode 
(green). The subplots present the transformation of the held 
prohles at the points corresponding to the vertical lines la¬ 
beled v-z indicated in Fig.[6l In each subplot the value of Eq 
{see Eq. ( [equation] [28] [] 28| in Ref. [18] for its dehnition} is 
identical for both modes. 


doubly degenerate. This means that to one value of the 
effective index (and the corresponding power density) 
correspond two solutions localized on one of the two core 
interfaces [compare green curves in Figs. I^a) and (e) for 
the ASl mode, Figs.j^a) and (e) for the AS2 mode and 
the green and the gray prohles in Fig. |^f)] 

The higher order modes with nodes resemble the modes 
of a linear slot waveguide with a high er refractive index 
than the one used here (see Sec. pg. Only symmetric 
(SI, S2, ...) and antisymmetric (ANl, AN2, ...) modes 
with nodes exist. Their dispersion curves are presented in 
Fig. [^and their held prohles are shown in Figs. ia)-(d). 
The dispersion curves of the modes with nodes start for 
/3 = 1 and their effective index grows with the increase 
of the power density Pc. 


In Fig. we present the dispersion relations obtained 
using the IM in a different coordinate frame. This time 
we use the total electric held intensity at x = 0 (the inter¬ 
face betw een the NSW core and the metal cladding) Eq 
{see Eq. ( [[equation] [28][]28 ) in Ref. HE]}. This quantity 
is one of the input parameters of the IM. The disper¬ 
sion diagrams P{Eq) have a drastically different charac¬ 
ter from the /3{Pc) diagrams presented in Fig. The 


difference is caused by the fact that Eq is a local quan¬ 
tity, whereas Pc is a global quantity, that results from 
the integration over the core width. In the coordinates 
of Eq^ the dispersion curves of the asymmetric modes are 
not degenerate. In Fig. we notice that for asymmet¬ 
ric modes, to a given value of (3 correspond two values 
of Eq^ that represent solutions localized on the left and 
right interface of the waveguide core. 

In Fig. the comparison of the field profiles of the 
three main modes is presented during their transforma¬ 
tion along the dispersion curve associated to the increase 
of Eq. The field profiles of the SO and ANO modes do not 
change qualitatively. On the contrary, the field profile 
of the ASl mode undergoes a qualitative transformation. 
For low Eq values, this mode is highly asymmetric and 
strongly localized on the right core interface x = d. With 
the increase of Eq the asymmetric profile becomes more 
symmetric, and at the point of bifurcation it perfectly 
overlaps with the symmetric mode [see Fig. |^c)]. For 
Eq values above the bifurcation point, the mode becomes 
asymmetric and it tends to localize on the left interface. 

In Fig. [^a similar transformation is shown for the SI 
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and AS2 modes. Here, with the increase of the peak 
amplitude of the soliton i^peak = Hy{x = d/2) first de¬ 
creases (it is the lowest at the bifurcation) and then in¬ 
creases, while side lobe peak amplitude i^iobe (located 
at X = 0 and x = d) oi the symmetric mode increases 
monotonously with Eq. The ratio i^peak/H^iobe first de¬ 
creases (up to the bifurcation) and then increases. In 
case of the asymmetric mode AS2, with the increase of 
Eq the soliton peak shifts from left to right. At the same 
time the amplitude of the left (right) side lobe increases 
(decreases). More results on the mode transformation, 
including the results obtained using the JEM are pre¬ 
sented in Ref. [20] . 
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FIG. 5. Typical magnetic field profiles Hy{x)^ obtained using 
the IM, corresponding to different dispersion curves indicated 
in Fig.[^ Abbreviations next to the subfigure labels indicate 
the dispersion curve to which a given profile corresponds. The 
color of the profile allows us to distinguish the mode symme¬ 
try: symmetric (S — blue), antisymmetric (AN — red), and 
asymmetric (AS — green). For asymmetric doubly degener¬ 
ate mode ASS the second profile is shown in gray. 


B. Single interface limit 


In Sec. 


section 


4][]IV in Ref. [18], describing the theo¬ 


retical derivation of the models for NSWs, we mentioned 
that in the limiting case, w here the integration con stants 
Co in Eqs. ( equation 8 ^ and ([equation 10 Toj ) or Co 


in Eqs. ( [equation][26] []26| ) and ([equation] [33][]33) (all 
equation numbers correspond to Ref. W are equal to 
zero, we recover the case of a single interface between 
a metal and a nonlinear dielectric. Looking at the field 
profiles of highly asymmetric modes ASl [see Eigs. [^a), 
(e)], we see that these modes are mostly localized at one 
interface only. Therefore, they can be well approximated 
by a solution of the single-interface problem. 


In Eig. we present the dispersion curves for the 
NSW obtain ed using the JEM [^{Ho)] {see Eq. ( |[eqna-| 
|tion] [1] [9]9a| ) in Ref. m} and the IM [/3(£’o)] (com- 
pare with Eig. [^. Additionally to the antisymmetric 
(red), symmetric (blue), and asymmetric (green) dis¬ 
persion curves, black dispersion curves obtained using 
single-interface models are presented. In the case of the 
JEM, the single-interface approximation was obtained 
using the ’field based model’ for configurations with semi¬ 
infinite nonlinear medium described in Ref. [21]. This 
model was used for a single interface between a metal 
and a nonlinear dielectric with the same parameters as 
our NSW. In case of the IM, the corresponding single¬ 
interface approximation was obtained using the ’exact 
model’ for configurations with a semi-infinite nonlinear 
medium described in Ref. [21]. 


In Eig. we see that both for the JEM and for the 
IM, the single-interface dispersion curve always lays be¬ 
tween the antisymmetric ANO curve and the symmetric 
SO curve. Eor high values of Eq, the asymmetric ASl 
curve becomes very close to the black curve, but remains 
slightly above it. The fact that the black curves over¬ 
lap with the green ASl curves confirms that the highly 
asymmetric ASl modes (for high effective index /3) are 
well approximated using the single-interface approach. 

Now, instead of using the corresponding models for 
configurations with a semi-infinite non linear medium, we 
will use the formulas found in Sec. [[section] [4] []IV 


Ref. [18] that give us the analytical expressions for the 
dispersion relations for the single-interface problem. In 
the case of the JEM, the analytical formula for the dis¬ 
persion relation of a single metal/nonlinear dielectric in¬ 
terface problem is given by Eq. ([equation] [42] [] 42) in 
Ref. [18]. In this equation, as in the entire formula¬ 
tion of the JEM, the primary parameter is the magnetic 
field amplitude at the interface Hq. Therefore, we are 
able to show the dependency described by Eq. (|[equa-' 


jtion] [42] []42| ) only in the coordinates where the effective 
index is presented as a function of the magnetic field am¬ 
plitude at the interface Hq [see yellow line in Eig. [^a)]. 
We observe that the dispersion relations calculated us¬ 
ing the field based mod el (balck curve) an d the yellow 
curve described by Eq. ( [[equation] [42][ 4^ overlap per¬ 
fectly. The single-interface dispersion curve, which cor¬ 
responds to the limiting case cq = 0 divides the dispersion 
plot P{Ho) into the regions corresponding to the node¬ 
less family and the family w ith nodes as predicted in Sec¬ 
tion [[subsection] [2][3]III B[ in Ref. [18] Above the cq = 0 
curve (for negative values of the integration constant cq), 
only node-less solutions exist. Below the cq = 0 curve 
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FIG. 6. Dispersion diagram for the symmetric NSW obtained using the IM. The vertical black lines indicate to the values 

of the total electric held amplitude Eq corresponding to the held prohles depicted in Figs. and Open circles correspond to 
the held prohles shown in Fig. 




Eo [V/m] 


persion relation for the single-interface problem is given 
by Eq. ( [equation] [46] []46| ) in Ref. [18]. The effective in¬ 
dex of the mode expressed as a function of the material 
parameters of the structure and the total electric held 
inte nsity at the interfa ce Eq. The curve described by 
Eq. ( [[equation] [46] []46| ) is plotted in yellow in Eig. [^b) 
and it overlaps well with the black curve obtained using 
the exact model. 

In case of the IM, the numerical results also show that 
the dispersion curves are divided in two families: with 
nodes and node-less. The regions of the dispersion dia¬ 
gram corresponding to these two families are separated 
by the curve described by the equation Co = 0 (black 
and yellow curves for the single-interface problem). In 
the frame of the IM, we could not prove this property 
analytically because the held plots in the IM are calcu¬ 
lated numerically. 

In Eig.[^ in the region of high ehFective indices, the dis¬ 
persion curves of the ASl mode overlap with the curves 
obtained using single-interface approximations. This 
conhrms our hypothesis that highly asymmetric modes 
ASl can be approximated by solutions obtained using 
the corresponding single-interface models. 


FIG. 7. Dispersion diagrams (a) /3(i7o) obtained using JEM 
and (b) /3(Eo) obtained using the IM for the symmetric 
NSW. Dispersion curves presenting single-interface approx¬ 
imation obtained using models derived in Ref. m are shown 
in black. Additionally, the curves corresponding to the an¬ 
alytical expression for the single-interface dispersion relation 
{Eqs. ( [[equatioiT [42] []4^ and ( [equation] [46] []4^ in Ref. [TS] } 
are shown in yellow. 


(for Co >0), only solutions with nodes exist. 

In case of the IM, the analytical formula for the dis- 


C. Comparison with linear states 

In Sec. m while discussing the field profiles of the 
modes belonging to the family with nodes, we noticed 
that they resemble higher-order modes of linear slot 
waveguides with parameters similar to the NSW stud¬ 
ied here. In this section, we will explain the origin of the 
similarities between these nonlinear and linear modes. 

In Eig. we present the nonlinear dispersion diagram 
obtained using the IM for our NSW. In this plot, the ef¬ 
fective index of the mode [3 is presented as a function of 
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the averaged nonlinear index modification in the waveg¬ 
uide core (An): 

(An) — ^ J ^ J ( 2 ) 

where the nonlinear parameter n^ 

In addition to this plot we present also a dispersion 
relation (black curves in Fig.[^ of a linear slot waveguide 
with a homogeneous and linear core and the following 
parameters: ei = 63 = —90, n = no + Anun = 3.46 + 
Aniin and d = 400 nm. The parameters ei, 63 and no = 
y/ei ^2 are identical to these in the case of the nonlinear 
waveguide studied here. 

We notice that this linear dispersion diagram is simi¬ 
lar to the dispersion plot of the NSW. For the core with 
index n = no only two modes are present and they are 
the linear counterparts of the modes SO and ANO. With 
the increase of the core index n, the effective index of 
these modes increases and they become closer to each 
other. At Anun ~ 0.1, a higher-order linear mode ap¬ 
pears that is a counterpart of the SI mode. For Anun ~ 2 
and Aniin ~ 3.5, another two higher-order modes ap¬ 
pear. They are the linear counterparts of the ANl and S 2 
modes, respectively. The effective index of these modes 
increases rapidly with the increase of Anun- The only 
modes not present in the linear dispersion curves are 
the asymmetric modes ASl, AS 2 , ...and the symmet¬ 
ric node-less modes SI, SII, etc. The asymmetric modes 
can not be observed in the linear case because nothing 
breaks the symmetry in the symmetric linear slot waveg¬ 
uide. The node-less symmetric modes are not supported 
by the linear slot waveguide because they have purely 
nonlinear solitonic character [see Figs. and [^e)]. 

The dispersion curves of the nonlinear modes ANO 
and SO overlap with the corresponding linear dispersion 
curves only for small (An) values. The nonlinear modes 
increase their effective indices (3 faster than the linear 



0.001 0.01 0.1 1 4 


<An> 

FIG. 8. A comparison of the nonlinear (red, blue and green 
curves) and the linear dispersion plots (black curves) of the 
symmetric slot waveguides. In case of the linear waveguide 
(An) is equivalent to Anun. Points correspond to the modes 
presented in Fig. 
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FIG. 9. A comparison of (a), (c) Hy{x) and (b), (d) Ez{x) 
for the nonlinear modes SI (blue curve) and ANl (red curve) 
and the normalized profiles of their linear counterparts (black 
curves) at the common points of the black (for linear slot 
waveguide) and blue or red (for NSW) dispersion curves in¬ 
dicated by open circles in Fig. 


modes. In case of higher-order modes SI, ANl and S 2 , 
the dispersion curves of the linear modes lay below the 
corresponding nonlinear modes. There is only one com¬ 
mon point per mode for these curves (indicated by an 
open circle in Fig. and it turns out that at this point, 
the index distribution induced by the nonlinear mode in 
the nonlinear core is flat (data not shown). 

Figure present the comparison of the field profiles 
Hy (x) and Ez{x) for nonlinear SI and ANl modes and 
their linear counterparts, at the points where the index 
distribution induced by the nonlinear mode in the non¬ 
linear core is fiat. We observe that the nonlinear profiles 
overlap perfectly with the profiles of the linear modes nor¬ 
malized to the same amplitude as the nonlinear modes. 

The results presented here prove that the modes with 
nodes found in the NSW are close to the modes of the 
linear slot waveguide with similar opto-geometric param¬ 
eters. We explain the similarities between these nonlin¬ 
ear and linear modes using the self-coherent definition 
of nonlinear modes. This definition was introduced by 
Townes and co-workers in Ref. [ 22 ] and was used later in 
other works (e.^.. Ref. [23|). It defines a nonlinear mode 
as a linear mode of a linear (graded refractive index) 
waveguide that is induced by the light distribution of this 
mode. According to this definition, there is no difference 
between the nonlinear modes of the NSW for which the 
nonlinear index modification has a flat distribution and 
the linear modes of the waveguide with higher, uniformly 
distributed refractive index of the linear core. 
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FIG. 10. Average nonlinear index change at the appearance 
of the asymmetric ASl modes (An)th as a function of the 
absolute value of (a) the metal cladding permittivity of the 
symmetric waveguide |ei| = jeal and (b) the linear part of the 
nonlinear core permittivity ei ^2 All the other parameters of 
the NSW are identical to these used in Sec. ura 


D. Permittivity contrast 


In Ref. [14] , we have studied the influence of the width 
of the NSW core on the nonlinear dispersion for this 
this structure. Here we will discuss the influence of the 
permittivity contrast between the dielectric core and the 
metal cladding on the nonlinear dispersion diagrams of 
symmetric NSW. 

First, we will discuss the influence of the metal 
cladding permittivity on the nonlinear dispersion dia¬ 
grams of NSW. We have studied the dispersion plots 
for the NSWs with identical parameters as these used in 
Sec. |II A| but with different values of the metal cladding 
permittivity. We observe that the cladding with higher 
permittivity (lower in absolute value) allows us to reduce 
the (An) threshold values where the bifurcation of the 
ASl mode occurs. For metals with permittivity equal 
to —40, the bifurcation occurs at (An) « 0 . 02 , which is 4 
times lower than in the case of ei = 63 = — 90. For met¬ 
als with permittivity equal to — 20 , the bifurcation occurs 
at (An) ^ 8 • 10 “^, which corresponds to the reduction 
of the bifurcation threshold by two orders of magnitude 
with respect to the configuration with ei = 63 = —90. 
The dependency of the ASl mode bifurcation threshold 
(An)th on the metal cladding permittivity is illustrated 
in Fig.[lQ|^a). Looking at this plot, we conclude that with 
the increase of the metal cladding permittivity (decrease 
of its absolute value) the bifurcation threshold of the ASl 
mode decreases. This decrease is slow in the range of high 
index contrast between the metal and the nonlinear di¬ 
electric permittivity. Although, for |ei| = |e 3 | close to € 1^2 
the decrease of the bifurcation threshold is more rapid. 
Changing the metal permittivity from —20 to —15 al¬ 
lows us to decrease the bifurcation threshold by almost 
two order of magnitudes. For the metal cladding permit¬ 
tivity ei = 63 = —15 the bifurcation threshold is at the 
level of (An) ^ 10“^. This is four orders of magnitude 
lower than for the ei = 63 from the range [-400,-90], 
for which the bifurcation occurs at (An) ^ 0 . 1 . 

Next, we will study the influence of the change of the 


core permittivity on the dispersion diagram of the sym¬ 
metric NSW. We analyzed the plots of the dispersion 
curves for the NSWs with different linear parts of the 
core permittivity e/^ 2 - All the other parameters are iden¬ 
tical to these used in Sec. EU The behavior of the bi¬ 
furcation threshold expressed as the averaged nonlinear 
index modification (An) is presented in Fig, [l^b). The 
increase of the linear part of the core permittivity € 1^2 is 
accompanied by a monotonous decrease of the bifurcation 
threshold. From Fig. we notice that the increase of 

€ 1^2 from 1 to 25 results in the decrease of the bifurcation 
threshold by approximately three orders of magnitude. 

It is interesting to remind that, in the case of changing 
the permittivity contrast by varying the metal cladding 
permittivity [see Fig. [TQ|[a)], we observed a decrease of 
the bifurcation threshold for the ASl mode with the de¬ 
crease of the permittivity contrast between the cladding 
and the core permittivity. On the contrary, decreasing 
the permittivity contrast by changing the core permit- 
tivty, leads to the increase of the bifurcation threshold 
[see Fig.fT^b)]. 

This phenomenon can be explained using the field pro¬ 
files of the symmetric mode for different values of core 
and metal permittivities. We observe that increasing 
the permittivity of the core or increasing the permittiv¬ 
ity of the metal (decreasing its absolute value) leads to 
symmetric modes that are more localized on the waveg¬ 
uide interfaces and look closer to separate plasmon on 
both metal/dielectric interfaces. Because the overlap and 
therefore the interaction between the two plasmons is 
weaker, it is easier to break the symmetry of the mode. 
This explains the decrease of the bifurcation threshold. 
We conclude that changing the permittivity contrast by 
varying the linear part of the nonlinear core permittiv¬ 
ity, has opposite effect than changing the permittivity 
contrast by varying the metal cladding permittivity. 


III. RESULTS FOR ASYMMETRIC 
STRUCTURES 

In Sec. [TT| we have comprehensively discussed dis¬ 
persion diagrams and mode profiles in symmetric NSW 
structures. In this section, we will discuss the influence 
of the NSW asymmetry on the dispersion curves. The 
asymmetry is introduced by sandwiching the nonlinear 
core by metals with different values of the permittivity 
on both sides. Asymmetric NSW structures have not 
been studied before in literature. Here we present the 
analysis of these structures for the first time. 


A. Dispersion relations 


Figure 11 presents the nonlinear dispersion diagram 
obtained using the IM for the structure with the following 
parameters: core permittivity e /^2 = 3.46^, the second- 

order nonlinear refractive index = 2 • 10“^^ m^/W, 








core with d = 400 nm, metal permittivities ei = —110, 
63 = —90 at a free-space wavelength A = 1.55 jim. These 
parameters are identical to those for the structure studied 
in Sec. IIA except for the metal permittivities. Here the 
permittivity of the left metal layer is decreased to —110 
making the structure asymmetric. 



10 ^ 10 ® 10 ® 10 ^® 
Pc (W/m) 


FIG. 11. Dispersion diagram obtained using the IM for the 
asymmetric structure with ei = —110 and 63 = —90 (for the 
scheme of the structure see Fig. [^. Blue curves correspond 
to the modes for which sgn[F^x,o] = sgn[F^x,(i] and red curves 
correspond to the modes for which sgn[F^a,,o] = — sgn[F^a,,(i] 
{see Eq. ( | [equation][28][] 2^ in Ref. [TS] for the notations}. 
Compare this dispersion diagram for the asymmetric struc¬ 
ture with the dispersion diagram for the symmetric structure 
presented in Fig. [^b). 



FIG. 12 . Dispersion curves /3(Eo) for the asymmetric struc¬ 
ture with ei = —110 and 63 = —90. Compare this dispersion 
diagram for the asymmetric structure with the dispersion di¬ 
agram for the symmetric structure presented in Fig. 


In the asymmetric structure only asymmetric modes 
are present. However, in the dispersion diagram shown 
in Fig. we divide the modes in two groups: modes 
that resemble the antisymmetric modes of the symmetric 
structure for which sgn[F^a,^o] = “ ^^^[Ex,d] (red curves 
labeled AN-like) and modes that resemble the symmet¬ 
ric or asymmetric modes of the symmetric structure for 
which sg n[F^a;,o] = sgn[F^a;.d ] (blue curves labeled S-like) 
{see Eq. ([equation] [28] []28) in Ref. [18] for the notations 


of the electric field components}. 

We compare the nonlinear dispersion curves for the 
asymmetric structure presented in Fig. [E with the dis¬ 
persion curves obtained for the symmetric structure 
shown in Fig.[^ We notice that the dispersion curves for 
the symmetric and antisymmetric modes from the family 
with nodes did not change much. The number of modes 
and the character of their dispersion curves is conserved. 
The main difference between the dispersion curves of the 
asymmetric and symmetric structures can be observed 
for the symmetric and asymmetric modes of the node¬ 
less family. The asymmetry of the structure lifts the dou¬ 
ble degeneracy of the asymmetric branch ASl (see green 
curve in Fig.[^. This branch splits into two branches (see 
Fig. [IT] ). One of them (the branch with lower effective 
indices /3) is a continuation of the symmetric-like funda¬ 
mental mode (blue curve) that starts for small power den¬ 
sity Pc levels. The second branch lays along the first one 
but has slightly higher power levels (branch with higher P 
values). The degeneracy of the higher-order asymmetric 
modes is also lifted by the asymmetry of the structure. 
These branches also split into two separate branches, sim¬ 
ilar to the case of the ASl mode. It is difficult to observe 
this effect in Fig.[^ where the power density in the core 
is used as abscissa (even enlarging the region of interest), 
because the two dispersion curves into which the disper¬ 
sion curve of the higher-order asymmetric mode splits lay 
very close to each other. The degeneracy lift of the AS 2 
mode can be however observed from the dispersion curve 
I3{Eq) presented in Fig. 
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where the effective index is 
shown as a function of the field intensity at the left core 
interface. In these coordinates the separation of the SI 
and AS 2 curves reflects the degeneracy lift of the AS 2 
mode. 


B. Permittivity contrast study 


To finish our discussion of the asymmetric NSW prop¬ 
erties, we directly compare the dispersion diagrams ^{Pc) 
of the symmetric structure with these of the asymmetric 
structures. In Fig. the dispersion plot of the sym¬ 
metric structure (ei = 63 = —90, see Fig. is com¬ 
pared with the dispersion plot for the asymmetric struc¬ 
ture (ei = — 110,63 = —90, see Fig. 11). Only a vicinity 
of the bifurcation point of the ASl mode is presented. 
We observe that, for low Pc values, the dispersion curves 
of the two low-power modes are slightly modified due to 
the waveguide asymmetry. For higher values of Pc, the 
dispersion curve of the fundamental mode (upper blue 
curve) exactly overlaps with the dispersion curve of the 
asymmetric mode of the symmetric structure. This is 
a consequence of the fact that the field profiles corre¬ 
sponding to this upper blue curve are strongly localized 
on the interface with the metal with higher value of the 
permittivity. These profiles resemble the profiles of the 
highly asymmetric modes of the symmetric structure [see 
Fig-i a)]. Therefore, we are not surprised that these two 



















9 


dispersion curves overlap. The second curve that results 
from the degeneracy lift of the asymmetric mode lays 
above (in terms of Pc) the dispersion curve of the asym¬ 
metric mode ASl (green curve). 



Pc (W/m) 

FIG. 13. Dispersion curves of the asymmetric NSW with ei = 
— 110,63 = —90 (blue curves) and the symmetric structure 
ei = 63 = —90 (green curves). 



FIG. 14. Dispersion curves of the asymmetric NSWs with 
6 i = —70 and 63 = —90 (blue curves), 61 = —50 and 63 = —90 
(red curves), and the symmetric structure 61 = 63 = —90 
(green curves). 

In Fig. we present a comparison of the dispersion 
curves of the symmetric structure (ei = 63 = —90, see 
Fig. and the asymmetric structures, where one of the 
metal permittivity values is higher than in the case of 
the symmetric structure. The dispersion curves of the 
symmetric structure (green curves) are compared with 
these of the asymmetric structures with ei = —70, 63 = 
—90 (blue curves), and ei = —50, 63 = —90 (red curves). 

In the case illustrated in Fig. contrary to the one 
presented in Fig. it is the higher (in terms of Pc) of 
the two curves that result from the lift of the degeneracy 
that overlap with the dispersion curve of the asymmetric 
modes of the symmetric structure. This higher curve cor¬ 
responds to the modes that are localized on the interface 
between the core and the metal with permittivity equal 
to —90. For the structures studied in Fig. e = —90 


is the lowest cladding permittivity. For that reason, the 
dispersion curves corresponding to the mode localized on 
the interface with metal with lower permittivity, overlap 
with the dispersion curves of the symmetric structure. 

Another effect that can be observed in Fig. is that 
with the increase of the structure asymmetry |ei — 63 ! the 
separation of the two curves that appear as a result of the 
degeneracy lift, increases, as expected. In the limiting 
case 61 ^ 63 , these two curves merge into one doubly 
degenerate curve. 


IV. STABILITY OF THE MAIN SOLUTIONS 
FOR SYMMETRIC WAVEGUIDES 

In the previous sections, we have studied the stationary 
properties of plasmon-soliton waves using two different 
modal approaches. From both theoretical and practical 
points of view, the issue of the stability of these waves 
arises. In several works, the general problem of the sta¬ 
bility of nonlinear waves was studied [24H26] . Despite an 
enormous interest in the properties of nonlinear waves 
over the last decades, there in no universal condition on 
their stability [191123]. most of the cases, the stability 
must be studied numerically for each of the cases sepa¬ 
rately. Stability of nonlinear guided waves in fully dielec¬ 
tric structures was studied numerically in Refs. 

|33]. 

In structures made of metals and nonlinear dielectrics, 
due to the presence of media with negative permittiv¬ 
ity, the problem of stability of plasmon-solitons is dif¬ 
ficult to study even numerically. Only in Refs, 
the stability of plasmon-solitons was analyzed for the 
single nonlinear dielectric/metal interface case, using nu¬ 
merical algorithms (like finite-difference time-domain — 
FDTD |36l|37]). The propagation of light in plasmonic 
couplers was studied using Fourier methods based on 
mode decomposition in linear [38] and nonlinear [39] 
regimes. In this section, we study the stability of the 
plasmon-soliton waves in symmetric NSWs using two 
methodologies: (i) the topological criterion for funda¬ 
mental modes of nonlinear waveguides derived in Ref. [23] 
and (ii) two numerical full-vector methods (using COM- 
SOL [40] and nonlinear FDTD implemented in meep [4TJ 

m)- 


A. Theoretical arguments 

We use here the topological criterion presented in 
Ref. [23] that is based on the linear stability analysis [43] 
and the Vakhitov-Kolokolov criterion [44]. The stability 
criterion presented in Ref. m, is based only on the topol¬ 
ogy of the nonlinear dispersion curves and the stability 
of the modes can be read by analyzing ^{Itot) diagrams 
in which /tot = J(^)dx where I(x) is the intensity 
density. The validity of this approach was confirmed in 
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p I 


possibly / 
stable 


unstable 


tot 


FIG. 15. Rules of assigning the stability of the modes for two 
specihc cases extracted from Fig. 2 in Ref. [23]: (a) the fold 
bifurcation (open circle) and (b) the Hopf bifurcation (open 
square). Thick lines indicate a doubly degenerate branch, 
whereas thin lines indicate non-degenerate dispersion curves. 


multiple settings dealing with purely dielectric structures 

mMm 

First, we will recall the principle used to estimate 
the stability of nonlinear modes using the criterion from 
Ref. [23]. Then we will use it to analyze the stability of 
some of the plasmon-solitons found in NSWs. 

The stability criterion derived in Ref. [23] uses sev¬ 
eral assumptions. It provides stability for the funda¬ 
mental nonlinear modes in structures composed of ar¬ 
bitrary nonlinear material distributed nonuniformly in 
the transverse direction. The derivation of the stability 
criterion from Ref. [23] is obtained in the weak guiding 
approximation for which the electric field satisfies the 
scalar wave equation. In our study of the TM polarized 
waves, we consider the case in which it is the magnetic 
field component that satisfies the scalar wave equation 
in Ref. [H]]. We are fully aware 


[Eq. ([equation] [5 


of the fact that the metal/nonlinear dielectric structures 
studied here, in which plasmon-soliton waves propagate, 
do not fulfill the weak guiding approximation due to high 
permittivity contrast between the metal and the nonlin¬ 
ear dielectric. This means that interesting nonlinear ef¬ 
fects will occur for quite high nonlinear permittivity mod¬ 
ifications. In spite of this fact, we use here the criterion 
from Ref. [23] , because the dispersion diagrams obtained 
for our structures have similar character to the dispersion 
plots of the fully dielectric structures where the criterion 
is applicable and because, as it will be shown below, two 
different numerical propagation simulations of the full 
vector nonlinear problem confirm at least partially the 
theoretical predictions. 

In Fig.[^ the rules that will be required here to deter¬ 
mine the stability of the modes derived in Ref. m are 
schematically shown. Consider the dispersion relation 
presented in Fig.[^ It shows a zoom of a dispersion dia¬ 
gram, using /tot as variable, for a region that contains the 
dispersion curves of the main modes for the same struc¬ 
ture as the one presented in Fig.[^ The stability of modes 
changes only at the bifurcation points [23] . To determine 
the stability, first we have to identify all the bifurcation 


points on the dispersion diagram /3(/tot)- In Fig. 16 



ItotlW/m] 


FIG. 16. Zoom on the region of the dispersion diagram with 
the birth of the hrst order asymmetric mode. Bifurcation 
points are marked with an open circle for fold bifurcation and 
an open square for Hopf bifurcation. The numbers facilitat¬ 
ing the stability analysis are assigned to the sections of the 
dispersion curves according to the rules presented in Fig. 
Labels ’ps’ and ’u’ denote possibly stable and unstable modes, 
respectively. 


the bifurcation points are located at the points where 
intensity /tot has its local minima or maxima (point indi¬ 
cated by an open circle — so-called fold bifurcation [45] ) 
or where another branch appears [point indicated by an 
open square — so called Hopf bifurcation associated with 
the birth of a doubly degenerate branch (to a single point 
on this branch correspond two asymmetric field profiles)]. 
Modes appear from or disappear at the points of bifurca¬ 
tion. The next step is to label the sections between the 
bifurcation points with numbers. The numbers are as¬ 
signed in the following way. At first, we arbitrarily choose 
one section and label it with any number (in Fig. l^we 


labeled the low intensity section of the symmetric dis¬ 
persion curve with a number 0). The numbers of all the 
other sections of dispersion curves are assigned using the 


geometric rules given in Fig. 15 


Finally, after having numbered all the sections of the 
dispersion curves, we can read the stability of the modes 
directly from the P{Itot) dispersion curves. The topolog¬ 
ical stability criterion presented in Ref. [23] tells us that 
only the modes corresponding to the parts of the curves 
with the largest number are possibly stable. In Fig. 
only the modes labeled by 0 are possibly stable (ps). All 
the other modes are unstable (u). The stability of all the 
possibly stable modes can be specified at once, as soon as 
the stability of one of them is determined. The stability 
can be determined either using numerical methods or the¬ 
oretical arguments. The low-intensity section of the sym¬ 
metric branch in the linear limit corresponds to a linear 
plasmon in metal/insulator/metal (MIM) configurations, 
which is stable. Therefore, the solutions corresponding 
to this section of the nonlinear dispersion curves should 


be stable. This hypothesis will be confirmed in Sec. IV B 


The high intensity section of the symmetric branch 
(above the Hopf bifurcation) corresponds to unstable so¬ 
lutions. On the contrary, the section of the asymmetric 
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branch just above the Hopf bifurcation should correspond 
to stable solutions, because the stability properties of the 
sections with the same number are the same [23] . On the 
asymmetric branch (at /3 ~ 10) another bifurcation oc¬ 
curs (fold type bifurcation indicated by an open circle). 
The high effective index section of the asymmetric branch 
(above the fold bifurcation point) is unstable. 

B. Numerical simulations of nonlinear 
propagations 

In the previous section we provide some results about 
the stability of the plasmon-solitons of the lowest orders 
using the topological criterion derived in Ref. [23j. In 
the NSWs studied here the weak guiding approximation, 
used in the derivation of this topological criterion, is not 
fulfilled. This fact makes the conclusions drawn using 
the criterion not definitive. For this reason we also inves¬ 
tigate the stability by full-vector numerical simulations. 

First, we have used the capabilities of the FDTD 
method [31137] implemented in the meep software mi 
l42l [46] . The metal permittivity is described by a Drude 
model to obtain the fixed negative value used at the 
studied wavelength in the semi-analytical models we 
developped. The useful computational domain is sur¬ 
rounded at its four edges by absorber regions that pre¬ 
vent back-reflected fields more efficiently than the per¬ 
fectly matched layers that have also been tested during 
our FDTD simulations. An example of the FDTD prop¬ 
agation of the asymmetric plasmon-soliton is presented 
in Fig. [T^a) through the evolution of the electric field 
component Hy^ where the sinusoidal phase modifications 
are visible. 

This result provides a confirmation of the stability of 
the first asymmetric mode in the nonlinear slot waveg¬ 
uides when the intensity is above a critical threshold. It 
can be noticed in Fig. Hlfa) that the plasmon-soliton 
profile is not fully stationary. This is due to the fact 
that the used current excitations in our FDTD simula¬ 
tions do not generate perfectly the field profile of the 



2 [Mm] 


FIG. 17. Evolution of the Hy field profile of an asymmetric 
plasmon-soliton for the stable case (a) with (An) = 0.0138 
and for the unstable case (b) with (An) = 0.005. These simu¬ 
lations are realized using the FDTD method implemented in 
the MEEP software. The parameters are given in Fig. [T^ 


asymmetric plasmon-soliton. The used asymmetric ex¬ 
citations contain components that weakly excite the an¬ 
tisymmetric plasmon-soliton (which is studied later in 
this section). When the excitations match perfectly with 
the asymmetric mode profile the observed non-stationary 
behaviour disappears as shown later when the simula¬ 
tion results from the other used numerical method are 
described. The main symmetric plasmon-soliton is easier 
to excite in a simple way due to its symmetry property 
as it can be seen in Fig. [Ts] where a stable and stationary 
propagation is shown. 

In order to obtain a more general view on the stabil¬ 
ity of the main modes of the NSW in the frame of the 
FDTD method, we systematically studied the propaga¬ 
tion properties of the three main modes as a function 
of the spatially averaged refractive index variation (An) 
[see Eq. (§]: the first symmetric, asymmetric, and an¬ 
tisymmetric modes. Typically, three cases occur in the 
simulation results: 

• Case 1: The mode is visible during the entire sim¬ 
ulation duration. This case is, for example, the 
one encountered for the asymmetric mode above a 
given threshold (An) as it can be seen for example 
in Fig. [T^a) or in Fig. 

• Case 2: The studied mode is generated at the be¬ 
ginning of the temporal evolution but after some 
times it does not propagate anymore in a self¬ 
similar way. This case is the one encountered for 
the main asymmetric mode below a given threshold 
(An) as shown for example in Fig.[T7|[b) where only 
the most stable part of the propagation is shown. 

• Case 3: The investigated mode is not generated 
by the chosen current source (symmetric, antisym¬ 
metric or asymmetric) used to excite it, even at the 
beginning of the temporal evolution and in the sur¬ 
rounding of the source. This behaviour is observed 
for the asymmetric mode below the critical power 
or critical (An) associated to the Hopf bifurcation. 

It is one of the main advantages of the FDTD method to 
be able to simulate temporal evolution even in the case 
of unstable modes unlike the other method used later in 
this section. 

As it is shown in Fig. obtained from the FDTD 
simulations, we are able to build a dispersion diagram 


-0.5 - 

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 

FIG. 18. Evolution of the Hy field profile of a symmet¬ 
ric plasmon-soliton for the stable case with (An) = 0.0018. 
These simulations are realized using the EDTD method im¬ 
plemented in the meep software. The parameters are given 
in Eig. 


18 













12 


for the first symmetric and asymmetric modes taking into 
account their stability properties. The given [3 values for 
unstable modes, corresponding to the case 2 in the above 
paragraph, are the ones extracted from the simulations 
results in the stable initial part of the evolution. It is 
evident that for the case 3 above, no dispersion data are 
obtained. 

The stability properties of the asymmetric mode from 
the FDTD simulations differ from the ones deduced from 
the topological criterion given in the previous section for 
the stationary case. The asymmetric mode is not stable 
just above the bifurcation (see case 2 above) for some 
range of (An) (see the thin green curve in Fig. and 
then it becomes stable when (An) increases (see the thick 
green curve in Fig. 19). The instability of the asymmet¬ 
ric mode just after the bifurcation has already been de¬ 
scribed in the field of the spatial soliton studies mm- 
In our case, the instability can be observed in a relatively 
extended range of intensity or equivalently of (An). This 
extension of the instability could be due to the way the 
asymmetric mode is excited in our FDTD simulations 
and/or to the fact that the metal permittivity is disper¬ 
sive due to the used Drude model. 

It is worth noting that the FDTD dispersion curve for 
the asymmetric mode differs at high (An) from the one 
computed using the stationary IM: here the (3 values are 
smaller and the FDTD curve stays concave while the sta¬ 
tionary one is convex. Similar saturation effects in non¬ 
linear full-vector temporal simulations have already been 
described e.^., in Ref. m- From the FDTD implementa¬ 
tion we use, we can not conclude about the stability prop¬ 
erty at higher intensities than the ones shown in Fig. 



FIG. 19. Dispersion and stability results, obtained from 
FDTD simulations, for the first symmetric and asymmetric 
modes of the symmetric NSW. Thick curves denote a sta¬ 
ble propagation while thin ones denote an unstable propa¬ 
gation. The parameters are the following: core permittiv¬ 
ity ez ,2 = 3.46^, the second-order nonlinear refractive in¬ 
dex = 2-10“^^ m^/W, core thickness d = 500 nm, 
metal permittivities € 1=63 = —6 at a free-space wavelength 
A = I.O /xm. 


due to the limitations of the nonlinear treatment used 
(see reference [42]). Consequently, we can not check the 
stability properties around or above the fold bifurcation 


point described in Section IV A 


As it was expected from the previous section, the first 
symmetric mode is stable at low (An) or equivalently at 
low intensities (see the thick blue curve in Fig. [l^ . 

Its stability is lost slightly before the birth and the 
partial propagation of the asymmetric mode (see the thin 
blue curve in Fig. 19). For all (An) values tested above 
this transition region, the first symmetric mode is un¬ 
stable. It is worth mentionning that the stability of this 
symmetric mode is recovered numerically as soon as the 
symmetry is forced in the FDTD simulations prohibiting 
the appearance of asymmetric behaviour. 

The topological criterion given in given in Ref. [23| 
can not be applied to the first antisymmetric plasmon- 
soliton mode because it is valid only for fundamental 
modes. Therefore, the stability of this mode can only 
be infered from numerical simulations. The first anti¬ 
symmetric mode starts, in the low-intensity regime, from 
the stable linear antisymmetric plasmon and there is no 
bifurcations on its dispersion curve. Therefore, we ex¬ 
pect this mode to be stable. An example of this stable 
propagation is shown in Fig. [^ We observe no change 
of the field profiles during the propagation (the antisym¬ 
metric excitation used does not contain any symmetric 
component). 

In Fig. the dispersion curve of the antisymmet¬ 
ric plasmon-soliton is given together with its stability 
property. The antisymmetric mode is stable up to the 
maximum intensity that can be treated with the FDTD 
implementation we use. 

The stability properties of the three main plasmon- 
soliton modes in nonlinear slot waveguides are also veri¬ 
fied using the nonlinear propagation scheme implemented 
in the lastest version of the RF module of COMSOL Mul¬ 
tiphysics uni- This approach was successfully used to 
study the stability of solitons in lattices built of metals 
and nonlinear dielectrics [5QH5^ . This method is limited 
to the cases where the studied mode is stable since the 
iterative numerical method used to compute the fields do 
not converge in other cases. 

According to the conclusions drawn from Fig. in 
Sec. |IV A| and from the FDTD simulations, the low- 
power section of the symmetric branch corresponds to 


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 


z[Mm] 


FIG. 20. Evolution of the Hy field profile of a stable anti¬ 
symmetric plasmon-soliton with (An) = 0.0225. These sim¬ 
ulations are realized using the FDTD method implemented 
in the meep software. The parameters are the same as in 

Fig-lU 
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FIG. 21 . Dispersion and stability results for the first antisym¬ 
metric modes of the symmetric NSW obtained from FDTD 
simulations. The thick curve denotes a stable propagation. 
The parameters are the same as in Fig. |19| 



FIG. 22 . Evolution of the electric held norm during the 
propagation of the symmetric mode located below the Hopf 
bifurcation threshold. The average nonlinear index change 
in the core induced by this mode is equal to (An) = 10“^ 
and the propagation distance is approximately 13 free-space 
wavelengths. The parameters are the following: core per¬ 
mittivity 6^,2 = 3.46^, the second-order nonlinear refractive 
index n^^^ = 2 • 10“^^ m^/W, core thickness d = 400 nm, 
metal permittivities ei = es = — 20 at a free-space wave¬ 
length A = 1.55 /xm. These simulations are realized using the 
GOMSOL software. 


stable solutions. This result is also conhrmed by the 
simulation presented in Fig. obtained for a NSW 
with d = 400 nm. No stable symmetric solution was 
found above the (An) transition region using the method 
implemented in COMSOL conhrming the FDTD re¬ 
sults already obtained. The stability of the asymmetric 
branch above the bifurcation region observed in Fig. 
is confirmed by these numerical simulations as shown in 
Fig. Figure presents the evolution of the Ex elec¬ 
tric field component for the asymmetric solutions in such 
a case. 

Figure!^ shows the transverse profiles of the symmet¬ 
ric and asymmetric plasmons-solitons in the NSW. For 
each symmetry type, we compare the profiles obtained 
using the interface model (these profiles are used as input 



FIG. 23. Evolution of the electric field norm during the prop¬ 
agation of asymmetric modes located between the Hopf bifur¬ 
cation and the fold bifurcation. The average nonlinear index 
change in the core (An) induced by these modes is equal to 
(a) 2 • 10“^, (b) 3 • 10“^, and (c) 4 • 10“^. The shown prop¬ 
agation distance is approximately 12 free-space wavelengths. 
The parameters are the same as in Eig. |22[ These simulations 
are realized using the GOMSOL software. 
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EIG. 24. Evolution of the Ex held profile during the propa¬ 
gation of the solution presented in Eig. [23|b) for a slot with 
d = 400 nm. These simulations are realized using the GOM¬ 
SOL software. The shown propagation distance is approxi¬ 
mately 6 free-space wavelengths. 




EIG. 25. Gomparison of the |E| profiles obtained using the IM 
(and used as the input profiles in the GOMSOL based simula¬ 
tions) and cuts of the held evolution in the middle of the prop¬ 
agation range (z = 9 /xm — 6 free-space wavelengths) and at 
the end of the propagation (z = 18 /xm — 12 free-space wave¬ 
lengths) for (a) the symmetric nonlinear plasmon-soliton (see 
Eig. \22\ and (b) the asymmetric nonlinear plasmon-soliton 
[see W.[lib)]. 

in the COMSOL based propagation simulations) with the 
cuts of the prohles presented in Figs. [^and[2^b) . These 
comparisons validate the accuracy of the evolution simu- 
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lations and consequently the results obtained concerning 
the stability properties of the symmetric and asymmetric 
modes in the NSW. One can notice that the stationary 
behaviour is more clearly seen in the COMSOL based 
simulations than in the FDTD ones. This is due to the 
fact that in the former case we directly use as input the 
profiles provided by the interface model while in the lat¬ 
ter case we use excitation current sources to generate the 
fields that mimick the stationary field profiles. Since we 
are investigating nonlinear phenomena, it is not possible 
to use in the FDTD simulations a part of a linear waveg¬ 
uide to filter the needed profile as it is usually done in 
FDTD based linear studies [36] . 


V. CONCLUSIONS 

We have provided detailed results for the plasmon- 
soliton waves in planar slot waveguides with a finite- 
thicknees nonlinear dielectric core. In symmetric struc¬ 
tures, using the semi-analytical models we developped 
for stationary states, we haved investigated the proper¬ 
ties of the first main modes and reported new higher or¬ 
der modes including asymmetric ones that exist at high 
intensities only. We have also described complete dis¬ 
persion diagrams for these different modes as a function 
of various quantities including the total power, the field 
value at one interface between the metal and the non¬ 
linear core, and also the spatial average of the nonlinear 
refractive index change. We have proven that the total 
intensity or equivalently the spatially averaged nonlinear 
refractive index change corresponding to the Hopf bifur¬ 
cation threshold from the first symmetric mode to the 
first asymmetric mode can be reduced by several orders 
of magnitude with an increase of the permittivity of the 
core or of the metal cladding. We have also proven the 
versatility of our semi-analytical models studying asym¬ 
metric structures. For such structures, we have described 
the impact of the metal permittivity contrast that lifts 
the degeneracy of the doubly degenerated asymmetric 
mode providing a more complex dispersion diagram than 
the one of a symmetric structure. 

Concerning, the stability of the main symmetric and 
asymmetric modes, we have used an already derived 
topological criterion established only in the weak guid¬ 
ance approximation being fully aware that our structures 
lay beyond its validity range. Nonetheless, as shown by 


full-vector simulations, the topological criterion predicts 
correctly the principal stability properties of the main 
modes of the studied planar nonlinear slot waveguides. 
Using this criterion, we have shown that the asymmetric 
mode emerging through a Hopf bifurcation at a critical 
intensity is stable between this bifurcation and a fold bi¬ 
furcation located at higher intensity level. The stability 
of this asymmetric mode is lost at this fold bifurcation. 
On the contrary, the symmetric mode is unstable for all 
intensity levels above the Hopf bifurcation while it is sta¬ 
ble below. 

Using two different full-vector numerical propagation 
methods, we have studied the stability of the three main 
modes: the symmetric, asymmetric, and antisymmetric 
modes. We have shown that the asymmetric mode is 
stable above a critical intensity slightly larger than the 
threshold associated with the Hopf bifurcation computed 
for the stationary states from our semi-analytical models 
at least up to the maximum level of tested intensities. 
The symmetric mode is shown to be unstable slightly 
below and slightly above the Hopf bifurcation threshold, 
and to be stable at lower intensities. For all tested in¬ 
tensities, these results confirm qualitatively the results 
derived from the topological criterion even if quantita¬ 
tive differences exist. Finally, we have also proven nu¬ 
merically that the anti-symmetric mode is stable in the 
entire range of tested intensities. 

These stability results together with those about the 
decrease of the bifurcation threshold should facilitate the 
design of specific structures in order to make possible the 
experimental observation of these plasmon-soliton waves 
more than thirty years after their theoretical discovery. 
Future studies should be dedicated to the further reduc¬ 
tion of the bifurcation threshold and to the study of more 
sophisticated configurations. 
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